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Abstract 

We present a graphical analysis of the adiabatic connections underlying double-hybrid density- 
functional methods that employ second-order perturbation theory. Approximate adiabatic connec- 
tion formulae relevant to the construction of these functionals are derived and compared directly 
with those calculated using accurate ab initio methods. The discontinuous nature of the approx- 
imate adiabatic integrands is emphasized, the discontinuities occurring at interaction strengths 
which mark the transitions between regions that are: (i) described predominantly by second-order 
perturbation theory (ii) described by a mixture of density-functional and second-order perturbation 
theory contributions and (iii) described purely by density-functional theory. Numerical examples 
are presented for a selection of small molecular systems and van der Waals dimers. The impacts 
of commonly used approximations in each of the three sections of the adiabatic connection are dis- 
cussed along with possible routes for the development of improved double-hybrid methodologies. 

Keywords: Density-Functional Theory, Perturbation Theory, Hybrid Functionals, Adiabatic 
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I. INTRODUCTION 



Hybrid density-functional methodologies are now amongst the most commonly applied 
methods in quantum chemistry. The idea to hybridize density-functional theory with wave- 
function based approaches can be traced back to the work of Becke [l, 3- Motivated by 
the adiabatic connection (AC) formalism 3ru] Becke suggested that the introduction of a 
fraction of orbital dependent exchange could be beneficial and deliver improved performance 
over standard generalized gradient approximation (GGA) functionals. Becke initially con- 
structed a half-and-half functional containing 50% orbital dependent exchange and 50% of 
a density-functional approximation to exchange based on a simple linear model of the AC. 
Although the performance of this initial model was disappointing subsequent empirical opti- 
mization of the weight of orbital dependent exchange based on thermochemical performance 
showed that a weight in the range 20%-30% delivered improved performance over standard 
methods. Later, a number of theoretical rationales for weights of orbital dependent exchange 
in this range have been pu t forward js, 9], however, as recently noted by Cortona [h]] and 
utilized by Guido et. al 11] , similar arguments can be made for a number of different values. 

The widespread adoption of exchange-based hybrid functionals for molecular applications 
can be understood from their improved performance in applications such as thermochem- 
istry and the determination of equilibrium molecular structures. Perhaps more crucially, 
the extra computational effort associated with the evaluation of the orbital dependent ex- 
change is sufficiently modest so as not to hinder the application of hybrid approaches to 
a wide variety of chemical problems. Nonetheless, the introduction of orbital dependent 
exchange is not a panacea for all issues associated with standard GGA functionals. For 
many properties, or for molecules at geometries away from their equilibrium structure, the 
inclusion of orbital dependent exchange may be detrimental. For the simple example of 



the H 2 molecule see Ref. 



12]. Furthermore, it is well known that error cancellations play a 



crucial role when combining GGA exchange and correlation components 13] and the com- 
bination of 100% orbital dependent exchange with typical GGA correlation functionals is 
therefore ineffective. The construction of optimal hybrid exchange-based methods requires 
a tradeoff between improving the description of the exchange component whilst maintaining 
a reasonable description of correlation effects. 

A natural extension of hybrid approaches is to consider the hybridization of the correlation 
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energy in addition to exchange energy component, resulting in so-called double hybrids. In 
recent years a number of empirical approaches which consider the post SCF addition of 
second-order perturbation theory (PT2) terms have been developed 14j-|22l| . perhaps the 
most well known of which is the B2-PLYP functional of Grimme 14j . More rigorous routes 
to combine DFT and PT2 approaches have followed the suggestion by Savin 23J to employ 
long-range only wave function approaches in tandem with specially designed short-range 
density functionals. From a density-functional point of view this approach can be justified 
by use of a generalized AC as shown by Yang 241 ] . In this manner implementations have 
been constructed which combine most of the models of quantum chemistry with density- 

2(J, M0ller-Plesset (MP) 



functional theory including, configuration-interaction theory 
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perturbation theory [27|, coupled-cluster theory |28|, multi-configurational self-consistent 
field theory j^, 30] and N-electron valence state perturbation theory 31]. 

The advantage of range-separated approaches is that the division of labour between the 
wave function and density-functional methods can to some extent be controlled by the 
manner in which the two-electron integrals are modified, typically by using an error-function 



attenuation 



32 



33|. T 



lis has been shown, for example, to be effective in the treatment of 



dispersion interactions 



lis 


las 


27, 


34 1 



34] . The disadvantage of such an approach is the need to develop 



specialized complementary density functionals, although this task has been undertaken at 



the local density approximation (LDA) 23_|, |35| , GG A 
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371 ] and meta-GGA 



levels. 



Recently, Sharkas et al. [39[ have used a similar approach to consider the standard linear 
AC to hybridize DFT with PT2 for functionals in which the fraction a c of MP2 correlation 
energy equals the square of the fraction a x of Hartree-Fock (HF) exchange. Fromager {40 1 
has recently extended this approach to incorporate more flexible two-parameter double- 
hybrid energy expressions that satisfy the inequality a c < a^. Commonly used empirically 
optimized double-hybrid methods satisfy this inequality but do not in general have values 
of a c exactly equal to the a^. This approach allows the construction and analysis of en- 
ergy expressions, similar to those used in empirical double-hybrid approaches, on a more 
sound theoretical footing. In addition, insights into the relative values of the exchange and 
correlation weights may be obtained and give a rationale for the values typically obtained 
by empirical optimization [40]. One key advantage of considering a linear AC, in place of 
a generalized path, is that the complementary density-functionals required may be read- 
ily derived by the application of uniform coordinate scaling relations to existing standard 
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functionals |4jJ-|45|, see Refs. [39_|, |46( for examples of methodologies utilizing this idea. 

Whilst the (generalized) AC provides a legitimization for combining density-functional 
and wave function methodologies the efficacy of such approaches can only be measured by 
careful benchmarking. Given the plethora of possible combinations of density-functional and 
wave function components, and the associated opportunities for error cancellations, assess- 
ment of double-hybrid forms is certainly non-trivial. Since each double-hybrid functional 
expression may be thought of as a model for the AC in this work we take an alternative 
approach and investigate the quality of the underlying AC integrands. Recently, it has 



become possible to calcu 
initio methodologies 47- 



ate the ACs for atomic and molecular systems using accurate ab 



52| . Here we will employ the method outlined in Refs. 



49 



which utilizes Lieb's definition of the universal density-functional [53], to generate bench- 
mark AC integrands for comparison with those generated by two-parameter double-hybrid 
approaches. 

In this work we analyze the linear AC integrands relevant to double-hybrid methodologies. 
We commence in Section [XT] by introducing the theory of the linear AC. In particular, we 
consider the novel division of the linear AC into segments and their description using density- 
functional perturbation theory, which is directly relevant to double-hybrid approaches. In 
Section UTTl we introduce the approximations necessary to practically compute these AC inte- 
grands. Firstly we give a brief overview of the method to determine these quantities using ab 
initio techniques, which serve as a benchmark in this work. This is followed by a description 
of the route used to calculate the AC integrands relevant to the Ai-B2-PLYP type function- 
als introduced in Ref. [401. The relationship between these two-parameter double- hybrid 
ACs and the standard B2-PLYP method is then discussed in this context. In Section HVl we 
present the calculated ACs for a number of small molecules and van der Waals dimers using 
these approaches. The results highlight the challenges faced in constructing double-hybrid 
functionals and the insights provided by our approach are outlined. In Section [V] we make 
some concluding remarks and discuss prospects for future work in light of our findings. 

II. THEORY 

In this section we introduce the linear AC and discuss how it may be partitioned into a 
number of segments. This partitioning of the adiabatic integrand enables the consideration of 
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approaches in which density-functional and wave-function based approaches for describing 
the integrand can be mixed. The segmented integrand is then expanded through second 
order in density-functional perturbation theory, providing a rigorous framework for one- 
and two-parameter double-hybrids. 



A. The linear adiabatic connection 



Let us consider a physical density n, which is associated with the ground state wave 
function, for the Schrodinger equation 

f + W ee + V)\^/) = E\^), (1) 



where V — j dr v(r) n(r) is the corresponding local potential operator, T is the kinetic 
energy operator and W ee is the two-electron repulsion operator. It is well known that 
formally exact expressions for the exchange and correlation energies associated with this 
density can be obtained when considering a fixed-density linear AC 0-7] between the non- 
interacting Kohn-Sham (KS) system and the fully-interacting physical system described by 
Eq. p. 

Introducing a parameter A to modulate the strength of the electronic interactions we may 
write the auxiliary partially interacting Schrodinger equation 

T + XW CC + V x ^j\^ x ) = S x \^ x ), (2) 

where the local potential operator V x = j dr v x (r) h(r) is constructed such that the density 
constraint 

nyx(r) = n(r), < A < 1, (3) 

is fulfilled. In the A = 1 limit, v A (r) and \1/ A should therefore reduce to the physical v(r) 
and respectively. On the other hand, for A = 0, the KS potential and determinant $ KS 
are recovered. Using the Hellmann-Feynman theorem 

(^ x \W ee \V x )+ / dr^-(r)n(r), (4) 



dA x 1 ee| J d\ 
the ground-state energy of the physical system, described by Eq. (TJJ), may be expressed as 
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an integral over the interaction strength on the interval [0, 1] 

f 1 df u 
E = / ^_du + £° 
Jo di/ 

= ($ KS |f + V |$ KS ) + E R [n) + / VV^Jn] du, (5) 

Jo 

where -EhM = 1/2 // ^( r i)^( r 2)/ r i2 dr 1 dr 2 is the Hartree energy. The exchange-correlation 
integrand W^Jn] can be further decomposed into exchange and correlation contributions as 
follows 

VCM = w>] + w c >], 

W x >] = ($ KS |# CC | $ KS ) - E n [n] = £ X HF [$ KS ], 

W c >] = (^|Wee|^) - ($ KS |# CC |$ KS ), (6) 

where _E^ F [$ KS ] denotes the HF exchange energy expression calculated for the KS determi- 
nant. 

The expressions of Eq. ([6]), whilst implicit functionals of the density, n, are expressed 
as explicit functionals of the partially- and non-interacting wave functions, \E r " and $ KS , 
respectively. Since in this work we wish to distinguish clearly between density-functional and 
wave-function based expressions for the exchange and correlation integrands we introduce 
the following notation for the exchange-correlation integrand in terms of explicit density 
functionals 

VCM = E x [n] + A"[n], (7) 

where E x [n] is the exact KS exchange density functional with a value equal to £' X IF [$ KS ] in 
Eq. ([H]). The correlation contribution to the integrand as an explicit density functional may 
be determined by considering the correlation energy of a partially-interacting system, 

E X c [n] = [ X W:[n]dv, (8) 



o 



which can be expressed in terms of the usual correlation density-functional (A = 1 limit) by 



means of a uniform coordinate scaling in the density 4lM44| . 

E^[n] = X 2 E c [n 1/x ], 



(9) 

ra 1/A (r) = (1/A) 3 n(r/A). 
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Differentiating this expression according to Eq. ([8]) gives the well-known result 

An , , dEZ\n] 



2z/£ c [n 1/l ,] + z/ 2 ^— ^ (10) 



<9z/ 

for the correlation integrand as an explicit density functional. Again we note that although 
A£[n] and W^[n] are formally equivalent we use different notations to emphasize how these 
quantities are determined. This notational distinction will be helpful in further discussing 
the segmentation of the adiabatic integrand and the description of these segments by different 
theories. 



B. Segmentation of the adiabatic integrand 

Practical routes to determine the AC integrands usually proceed either by determining 
wave functions that fulfil the density constraint in Eq. and then evaluating the expressions 
in Eq. (jBJ) or by evaluating the expressions of Eqs. (J7J) and (jTUJ) using an approximate 
exchange-correlation functional. The former route can be used to provide ab initio estimates 
of the AC integrand, whilst the latter is appropriate to provide estimates corresponding to 
purely density-functional approaches. 

Since the two-parameter double hybrids examined in this work involve both wave-function 
and density-functional contributions, we examine a segmentation of the linear AC to combine 
both approaches. To achieve this we decompose the exact exchange integrand in two parts, 

W>] = £ X HF [$ KS ] x X [0M (u) + EM x T [X2>1] {v), (11) 

and segment the correlation integrand into three parts, 

w c >] = [em-eT^+KH-KM 

+ (W\W ee \V') -E nx [n*»}) xX [0 ,A l[ (^) 

+ (E x [n] - £ X HF [$ KS ] + AcM " 

+ (^ Al |Wy^ Al ) -£ Hx [n*Ai]) xX [XlM (iy) 

+A:HxI m (,), (12) 
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where < Ai < A2 < 1 and EhxM — En[n] + E x [n\. The indicator function is defined as 
follows 

f 1 v G A 

X A {v) = { . (13) 

[0 v i A 

As long as the AC is exact, meaning that the density constraint is fulfilled, the expressions 
for the correlation integrand in Eqs. (EJ), ( fTUl) and (I12j) are equivalent. The advantage 
of not simplifying further the latter in the first and second segments, [0,Ai[ and [Ai,A2[ 
respectively, lies in the fact that approximate descriptions of the AC involving both wave 
functions and density-functionals can then be formulated and calculated easily by relaxing 
the density constraint. Once this constraint is relaxed the first two pairs of terms in the 
parentheses relating to the first and second intervals will not compensate anymore, since 
n 7^ UiSfu. However, as we shall see, their contributions may be expected to be small under 
certain conditions. 

The first segment of the correlation AC in the interval [0, Ai[ has been constructed from 
the pure density-functional exchange-correlation expression in Eq. ([7]), where (i) The HF 
exchange based on the KS determinant has been substituted for the exchange density- 
functional energy according to Eq. (ITT]) and (ii) The correlation density-functional integrand 
Ap[n] is removed by subtracting A^n^] based on the wave-function density and then re- 
placed, by (\E f!/ |H / e e |\I/ I/ ) — £?Hx[^*"]- Even when the density constraint is relaxed the latter 
term is the dominant contribution to the correlation integrand in this segment and is purely 
wave function-dependent. In this respect, the first segment is predominantly a wave function 
one. In the second segment for the interval [Ai, A2[, the correlation integrand varies as A"[n]. 
The additional wave function terms simply ensure a continuous transition from the first to 
the second segment, which can be referred to as hybrid wave function/density-functional 
segment. The final segment [A2, 1] is described within density-functional theory. As shown 
in Appendix [A], this partitioning of the integrand leads, by integration over [0,1], to the 
exact energy expression initially proposed by Fromager [40] : 

E = (^\f + \ 2 W ee + V\^) + E^ 2 [n^A, (14) 
where the complementary density-functional contribution equals 

Et" 2 [n] = ^hxcM + (Ai - A 2 ) (EbM + A^[n]), (15) 



with 



^HxcM = (1 - Xi)EuM + E c [n] - E^[n) 



(16) 



and the wave function \I/ Al is obtained self-consistently as follows 

<- nnn{<tf|f + X 1 W ee + V\^) + E^ KC [n^y (17) 

C. Density-functional perturbation theory 

A perturbation expansion of the exact ground-state energy expressed in Eq. f lT^|) can 



be obtained when solving Eq. (I17j) within MP perturbation theory |39l . |40| . By analogy to 
the HF approximation, the zeroth-order wave function $ Al is obtained when restricting the 
minimization in Eq. f fTTj) to single determinant wave functions 



$ Al <- mm|($|f + \ 1 W CC + V\<5>) +Eh xc K]}. (18) 
Let us introduce a perturbation strength a and the auxiliary energy 



(19) 





w 

ee 











+a(A 2 -Ai) /lT>a , AlhTfa , AlX +£ Hx c [n*-.Ai], 



with 



^ a ' Al mm|(^|f + A!?7 HP [<l> Al ] + aAxW Al 
+£h 1 xc M}=£ q ' Ai , 



(20) 



where t r HF [$ Al ] is the HF potential operator calculated with $ Al , and the perturbation 
operator 

AiVV Al = \i(w ee - U HF [& Xl ]^), (21) 



is the scaled fluctuation potential [39]. It is clear, from Eqs. ( fTB"l) and (1201) . that in the a = 
limit, \l> a ' Al reduces to $ Al , while, according to Eqs. (|T4j) . (|T7|) and (|T9|) . the auxiliary energy 
becomes, for a = 1, the exact ground-state energy since \l/ a ' Al reduces to \I/ Al . We note that 



the perturbation theory presented in this work differs from the one of Sharkas et al. 39] by 
the auxiliary energy which, in their approach, reduces to E a,Xl . Its perturbation expansion 
through second order, 

E a M = E <f>)Xi + aE WM + a 2 XlE { ^ Xl + 0(a 3 ), (22) 

where 

E (o)\i + E V)>* = (^\f + X l W ee + V\^ Xl )+E^ xc [n^ 1 ], (23) 

will be used in the following. The perturbation expansion of the wave function, which is 
deduced from Eq. (1201) . is however the same. The derivation presented here complements the 



work of Fromager 



40J, where perturbation theory was formulated in terms of an optimized 



effective potential (OEP) instead of a density-functional one. Finally, note the normalization 
factor in front of the two-electron interaction expectation value in Eq. (Tl9~l) . which must be 
introduced since the intermediate normalization condition 



= 1, < a < 1, (24) 



54j-|56| that, in this case, the wave function can be expanded 



will be used. It has been shown 
through second order as follows 

|^«,Ai) = |$^) + a Ai|*lS Al ) + a 2 |* (2)Al )+0(a 3 ), (25) 

where \I/mp Ai is the analog of the usual MP1 wave function correction. According to the Bril- 
louin theorem, the density remains unchanged through first order, leading to the following 
Taylor expansion, through second order, for the density: 

n^M (r) = n $ A! (r) + a 2 5n {2)Xl (r) + C(a 3 ), (26) 



so that self-consistency effects in Eq. (1201) do not contribute to the wave function through 



first order 



54j. Non-zero contributions actually appear through second order in the wave 



function [56j . From the wave function perturbation expansion in Eq. (1251) and the interme- 
diate normalization condition, we obtain the orthogonality condition ($ Al |\&Mp Ai ) = and, 
CIS db result, the following Taylor expansion 

A 

($ Al |# ec |$ Al ) 

(27) 

+2a\ 1 (^\W ee \^l Xl ) + 0(a 2 ), 





w 

' ' ee 








\$ai,Xi\ 
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where, according to Eq. (I2TI) . the first-order contribution can be rewritten as 

(^\Wj^Hl Xl ) = ($ Al |W Al |^ Al > 

= Eg? 1 , (28) 

since \I/mp Ai contains double excitations only. In addition, according to Eq. fl26|) . the com- 
plement density-functional Hxc energy difference can be expanded through second order 
as 

(^Hxc 2 - ^Hxc) = (^Hxf - E iLe) [ n <£ A i] 

- 2 / dr (w-Sy) l "- lfe(2,iI(r) (29) 

As a result, we finally obtain from Eqs. ( I19p . fl22|) . f l27|) and ( 129]) the following Taylor 
expansion for the auxiliary energy 

E *MM = E (P)\iM + aE ^lM + a 2 E m^M + ( a 3) j ( 30 ) 



where 



^,A 8 = ^(0)A l+ ^_ 1 ^ s ) [n#Ai]j 

^WAi^ = S (i)Ai + ( Aa _ a 1 )($ Ai |w / cc |$ Ai ), 



(31) 



The exact perturbation expansion of the energy through second order is then obtained in 
the a = 1 limit, which gives, according to Eq. ( 123|) . 

£ = (^\f + X 2 W ee + V\<l>^) + E^ 2 [n^ 1 ] 



+ (a? + 2A 1 (A 2 -A 1 ))eS Ai 

+ ... (32) 
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This energy expression is associated with the segmented correlation integrand in Eq. (|T2|) 
whose perturbation expansion through second order is obtained when decomposing the aux- 
iliary energy as follows 

E *MM = E a,o,o + / du+ / dzy 

Jo d " Jh d " 



($ KS \f + V\<$> KS )+E Ji [n} + f WM&v 

Jo 



+ / W* MM ' U du, (33) 



o 



ot,\\ ,A 2 ,f 
c 



where the exchange integrand is expressed as in Eq. (TTTj) and, in the a = 1 limit, W, 
reduces to the exact correlation integrand Wc[n]. The particular case Ai = A2 = corre- 
sponds to the exact KS theory: 



(34) 



As shown in Appendix [Bj the additional terms in the right-hand side of Eq. (|33|) (first 
line) introduce many-body perturbation theory corrections to the correlation integrand, in 
the first and second segments, which leads to the following perturbation expansion when 
considering the a = 1 limit: 

W c [n\ = (E x [n)-E™[$ KS ]+A»[n)-A»[n^} 

+ (^\W ee \<S>») + 2vE { 2 Xl - #H*[n*,]) x l [oM (u) 



+ (^>]-£n$ KS ] + A c >] 

-A^fn^] - J dr 5A A 75n(r)[n^ri^ Al (r) 
+ ($ Al |# ce |$ Al )+2A 1 ^ Al 
-En^n^} - J dr 5E Hx /Sn(v)[n^x 1 }5n (2)Xl (r) 

+ x X [\iM[W) 

+ ... (35) 

Note that, in the first segment, the MP2 contribution has been linearized for convenience. 
Nevertheless, after integration, the exact MP2 correlation energy is recovered (see Eq. (IB4I) ). 
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From the exact integrand expression in Eq. (jT2l) it is clear that, in both first and second 
segments, the wave function has been expanded in MP perturbation theory through second 
order. The third segment, which is the pure DFT part, is not modified by the perturbation 
theory treatment. 

Interestingly, the second-order correction to the density only appears in the second seg- 
ment. This is due to the fact that, in the particular case Ai = A2 (that is when the second 
segment disappears), the 2n + l rule is fulfilled (55), As a result, second-order corrections 
to the density (and therefore to the wave function) are absent from the second-order correc- 
tion to the energy. In the general case, where Ai 7^ A2, the second-order corrections to the 
density introduce a discontinuity at Ai in the correlation integrand. This would in principle 
disappear when expanding the integrand in the two first segments to infinite order, provided 
that the perturbation theory converges smoothly of course, which might not be the case in 
practice [561 ] . Finally we stress that, in the first two segments, all quantities calculated with 
& u [y > 0) correspond to correlation effects as defined in KS-DFT. In other words orbital 
relaxations which make Q u differ from $ KS contribute to the correlation energy, exactly like 
in second-order Gorling-Levy perturbation theory (GL2) [58 ]. 

III. COMPUTING THE AC 

In this section we introduce the methods used to calculate the AC integrands. The 
most accurate approach calculates the exchange and correlation integrands introduced in 
Sec. Ill AI by means of ab initio methods. The integrands calculated in this manner will be 
used to serve as a benchmark for more approximate approaches. To analyze practical double- 



hybrid methods based on the Ai-B2-PLYP methods, introduced in Ref. [40j, an approximate 
formulation of the second-order density-functional perturbation theory presented in Sec. Ill CI 
is considered. For comparison a similar approach is also applied to determine an AC for the 
conventional B2-PLYP functional. 
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A. ab initio estimates of the AC integrand 



To calcu 
Refs 



49 



ate accurate ab initio estimates of the AC we employ the methodology in 



-|52|. Following Lieb 53j, we write the auxiliary energy as 



£ x [v] = inf F x [n] + J w(r)n(r)dr , 
which gives the ground-state energy for the auxiliary Hamiltonian 



(36) 



H x [v] =f + XW ei 



dr v(r) n(r). 



(37) 



The universal density functional F x [n] can be expressed as a Legendre-Fenchel transform 
(convex conjugate) to the ground-state auxiliary energy £ A [t>], 



F x [n] = sup 



£ x [v) 



v (r)n.(r)dr 



(38) 



where the maximization is over a complete vector space of potentials. See Refs. |59M6l| for 
reviews of this approach to DFT. In the present work we employ ab initio approaches to 
calculate £ x [v] accurately and hence determine the functional F x [n] accurately. We note 
that even for approximate theories in finite basis sets where £ x [v] may not be guaranteed to 
be concave in v the functional of Eq. ( 13 8 p may still be constructed in a well defined manner, 
being conjugate to the concave envelope of £ x [v] at a given level of theory, which is denoted 
£ A [i>]. The concave envelope provides an upper bound, £ x [v] > £ x [v], with equality when 
£ x [v] is concave in v. In the limit of a full configuration- interaction treatment of correlation 
and a complete one-electron basis set the exact universal-density functional is recovered. 
For practical calculations we employ the algorithm proposed in Ref. 49( and implemented 



in Refs. 



5qL 



52] for arbitrary interaction strengths. The key aspect of this approach is to 



introduce an expansion of the potential 



«b(r) = ^cxt(r) + (1 - AK cf (r) + ^ b t9t(r) 



(39) 



which allows for the use of analytic derivatives in quasi-Newton approaches to perform 
the optimization of Eq. (I38p and determine the potential expansion coefficients {b t }. This 
opens up the possibility to perform calculations on molecular systems to complement earlier 
approaches applicable to atomic species 



47 



scheme detailed in Refs. 



48] . Here we use the second order optimization 



49J with a truncated singular value cutoff of 10 6 on the Hessian 



14 



and a gradient norm tolerance of 10 -6 . The Fermi-Amaldi potential is employed for v Te { 
and we use the same basis set {gt} for the potential expansion as is used for the molecular 
orbitals. To determine the densities n we use the Lagrangian method of Helgaker and 



J0rgensen 



62- 



651]. where required, to calculate the relaxed density matrices. For the H 2 
molecule we perform calculations at the FCI level to determine £ A [t>], n, and F A [n]. For 
all other species considered in this work we employ the coupled-cluster singles-doubles and 
perturbative triples [CCSD(T)] |66| method with all electrons correlated. All calculations 
are performed with a modified version of the DALTON2011 program 67|. 

To make the link to the AC we note that the Lieb functional of Eq ( 13"8j) is equivalent to 
the Levy-Lieb constrained search functional for canonical ensembles 53 1 



F x [n] 



min TyH 

7— ¥n 



M A [0]7 



X.n 



(40) 



where the minimization is over all density matrices with density n and 7 A ' n is the minimizing 
density matrix. The interacting functional F A [n] can be related to the non-interacting 
functional via 



F x [n] = F°[n} + [ ^Mdv. 
Jo 



du 

Identifying F°[n] with T s [n], evaluating -^F^n] by differentiating Eq. 
the Hellmann-Feynman theorem we obtain the usual AC expression 



F x [n}=T s [n}+ [ V^[n]di/, 
Jo 



where the AC integrand can be decomposed into 



W£ x >]dz/ = XE H [n] + \E x [n] + E x [n] 



(41) 

0]) and employing 
(42) 

(43) 



where £7hM is the classical Coulomb energy. The exchange energy component is given by 



EM = TiW ee f^ n - E n [n 
and the correlation energy of Eq. (JSJ) can be calculated using the correlation integrand 

u,n £,0,n\ 



(44) 



W c >] = TrW ee (l u ' n -f' n )- 



(45) 



Each of the energy components and their corresponding AC integrands can be calculated at 
each interaction strength following the optimization of Eq. ( |38|) . thereby mapping out the 
AC. 
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B. Ai-B2-PLYP double hybrid integrand 



We consider in this section an approximate formulation of the density-functional pertur- 
bation theory derived in Sec. Ill Cl where (i) the energy expansion is truncated at second order 
(ii) Becke exchange and LYP correlation density-functionals are used (iii) density scaling in 
the LYP correlation functional is ne glec ted (iv) second-order corrections to the density are 



neglected. The Ai-B2-PLYP energy 



40] is thus recovered from Eq. (1321) . 



E MM = ($^|T + V|$ Al ) +E B [n i x 1 ] + A 2 ££ F [$ Al ] 

+ (1 - A 2 )£> 4 a 1 ] + (l - Xi(2X 2 - Ai))^^] 
+A 1 (2A 2 -A 1 )E^ Al , (46) 

where, according to Eq. ffTSl) . the orbitals are computed as follows: 
<f> Al <- min|($|f + V\$) + E H [n^) + Ai^ F [$] 

+ (1 - Ai)^[n«] + (1 - A?)i£ w [n.] J. (47) 

The A!-B2-PLYP energy expression is formally identical to the conventional B2-PLYP one. 
The fractions of HF exchange a x and MP2 correlation energy a c can be identified as 



a x = A 2 Ai = a x - y/al~a c 

<— > , (48) 

a c = Ai(2A 2 - Ai) A 2 = a x 

as long as the condition a c < a x is fulfilled, which is usually the case in conventional one- 
and two-parameter double hybrids 4o|. In the spirit of Eq. ( 1331) . the Ai-B2-PLYP energy 
can be rewritten in terms of an exchange-correlation integrand, 



= ($ \T + V\$°) + E R ln^}+ [ n&'^dz/, (49) 

Jo 

where the KS-BLYP density n$ used as reference is recovered in the particular case Ai = 
A 2 = 0. The corresponding KS-BLYP energy can be expressed as 

£0,o = fioft + v\^)+E n [n^)+ f E ^[n %0 ]+2uE^[n^])du, (50) 
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where the uniform coordinate scaling in the LYP correlation integrand 

Af'ln] = 2vE™[n 1/v ] + ^ dE "> 1/A , (51) 

has been neglected. By analogy with Eq. ffTTj) . we define the Ai-B2-PLYP exchange inte- 
grand in terms of the HF and Becke exchange energies, both computed with the KS-BLYP 
determinant $°, as follows 

= £ X HF [^°] x T [0M (u) + £> 4o ] x T M (u). (52) 

The associated correlation integrand 

y^XiM,u = w^M," _ yyA^ (53) 

can then be deduced from Eqs. ( pf9l) and ( 1501) . like in Sec. Ill CI Since 

= E™[^] - E p [n^] - 2A 1 £ C LY >^ 1 ] 

+2A lJ E© A \ (54) 

we finally obtain the following expression for the A!-B2-PLYP correlation integrand 

Wt M ' v = {KM - E*[n*,\ + - E^°] 

+2uE { 2 Xl + 2u(E^ p [n ¥) ] - £ c LYP [r^])) x X [0M (u) 
+ ^a 1 ,a 2 ,a- + 2{y _ x x )E^[n^ x X [XlM (u) 
+2iyE^ YP [n^}xl [X2A] (u), (55) 

where X~[ means A — > X±. From Eq. (I55I) we note that the Ai-B2-PLYP correlation inte- 

A<Ai 

grand is continuous in Ax, even though approximate wave function and density-functionals 
are used. Qualitatively the behaviour of the AC can be understood by neglecting the vari- 
ation of all terms depending implicitly on v: 

+2<[n,„]xI [Al)A2| (,) 

+2vE? p [n^} xX [A2il] (z/). (56) 

In the first segment [0, Ai[, the slope of the Ai-B2-PLYP AC curve is dominated by the MP2 
correlation energy of the auxiliary Ai-interacting system. On the other hand, in the other 
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two segments, the conventional LYP correlation energy dominates the slope. Curvature 
could be introduced into the approximate ACs by considering higher-order MP terms and 
introducing density scaling effects. In this work we do not consider higher order perturbation 
theory energies, however, the effects of density scaling will be investigated further in Sec. IIVI 
For that purpose we define from Eq. (1351) a Ai density-scaled B2-PLYP (Ai-DS-B2-PLYP) 
correlation integrand: 

-(A C LYP >^] - 2vE™[n^\)) x I [0M {v) 

+ (v&* + Af p > 4 o] - Ar- Al K„]) x l [XlM (u) 

+An % ]xI [A2!l] (,), (57) 

where the density-scaled LYP correlation integrand is defined in Eq. (15"TT) . Note that the 
second-order corrections to the density have been neglected. As a result, the approximate 
Ai-DS-B2-PLYP correlation integrand remains continuous in Ai. Moreover, for simplicity, 
the orbitals used in Ai-DS-B2-PLYP and Ai-B2-PLYP schemes are the same, which means 
that density scaling has not been taken into account in the self-consistent calculation of the 
orbitals. 

The integrands of Eqs. (155]) and (1571) provide an approximate description of the AC for 



which the reference density is the KS-BLYP one. According to Eq. (T4T1) and Ref. 40], 
the local potential v u , which ensures that the density constraint on the v- interacting wave 
function is fulfilled, is approximated here by 

v u (r) -> v(r) + (l-u)[ + xl f + (1 - v 1 ) \ 1 f ■ (58) 

\ on{r) dn{r) I on(r) 

Due to the Brillouin theorem the zeroth-order density remains unchanged through first 
order and, as illustrated by the Ai-B2-PLYP AC curves in Sec. ITV] it does not vary signifi- 
cantly on the segment [0, Ai[. 
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C. Conventional B2-PLYP double hybrid AC integrand 



In conventional B2-PLYP calculations, the energy is calculated as follows 

E ax ' ac = ($ ax ' ac |f + V\T*' ac ) + E n [n^c] + a x E» F [T x ' ac ] 
+ (l - a x )£>^c] + (1 - a c )E c LYP [n^x, ac ] 

-=(2)a x ,a c , . 

+a c E MP . (59) 

This expression is formally identical to the Ai-B2-PLYP energy, provided that the condition 
a c < a x is fulfilled, and so the mapping in Eq. (JJHJ) between (a x , a c ) and (Ai, A 2 ) exists. The 
difference between the numerical values of the B2-PLYP and Ai-B2-PLYP energies lies in 
the computation of the orbitals used, where the two parameters a x and a c , instead of one, 
like in the Ai-B2-PLYP scheme (see Eq. (I47p ). are involved: 

T x ' ac <- mm^\f + V\^) + E u [n^} + a x E^ F [^] 

+ (1 - a x )£ x B K] + (1 - a c )E c LYP K]|. (60) 

As a consequence, there is an ambiguity in the way the correlation integrand should be 
defined for B2-PLYP. Indeed, the B2-PLYP energy expression cannot be derived rigorously, 
from the density-functional perturbation theory in Sec. Ill C\ as long as the orbitals are 



calculated according to Eq. ( 160]) . In this case, the Brillouin theorem cannot be applied [401 ]. 
which is an important difference with the Ai-B2-PLYP scheme, so that single excitation 
contributions to the double hybrid energy should in principle be considered. Nevertheless, 
it is interesting for analysis purposes to construct analytically a B2-PLYP integrand that 
can be compared to the Ai-B2-PLYP one. The simple segmentation 



pOx,«c -prO,0 f a * d — v,u 2 f ax d — a x ,j/ (ci\ 

E = E + / —E dv - —E du, (61) 

Jo dz/ Ja c dz/ 

could be used, but then the connection to Ai-B2-PLYP would be lost, simply because differ- 
ent intervals are used. In fact, segmenting the B2-PLYP energy in the same manner as the 
Ai-B2-PLYP one is not trivial. The main reason is that, in the Ai-B2-PLYP scheme, Ai and 
A2 are independent parameters but a c and a x are not, since the former depends on both Ai 
and A2. On the other hand, in B2-PLYP, a c and a x are independent parameters. Considering 
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that, in the particular case a c = a\ or equivalently Ai = A2, B2-PLYP and Ai-B2-PLYP are 
identical (E ' = E u ' u ), we propose the following segmentation, by analogy with Eq. fH9|) . 



E^ C = ET' V + I —E di> 

dz/ 



+ r A (T y - e"^-^ 2 ) diA (62 ) 

y^-v^ dz/V 7 

^ ^2 

Note that the derivative di? ' /dz/ is integrated up to a x , which ensures that the orbitals 
are calculated with the fraction a x of HF exchange. This is an important difference with 
Ai-B2-PLYP for which this fraction equals Ai = a x — v/a| — a c instead. As a result, the 
corresponding fraction of MP2 correlation energy must be reduced from to a c , which is 
exactly what the third term in Eq. f )62|) is devoted to. Finally, since the particular case 
o x = etc = corresponds to a standard BLYP calculation, the conventional B2-PLYP energy 
can be rewritten, according Eqs. fl50l) and fl62|) , and Appendix O, in terms of an exchange- 
correlation integrand 

E a *' ac = (&\f + V\$°) + E ll [n ¥) }+ ['W^du, (63) 

Jo 

where the exchange part is defined, like in the Ai-B2-PLYP scheme, from the HF and Becke 
exchange energies calculated for the KS-BLYP determinant 

Wr = £ X HF [^ ] x l [0M (u) + £>* ] x Z [fcil] („), (64) 
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and the associated correlation integrand equals 



xX r- — hi) 

+2(a x - \/a\- a^E^'""" 
+2,(CKo]-C[^l) 
+2(a x -u)E^ p [n_ a ^_ {a ^ )2 ] 

xX r 1 h>) 

+2,Ck.]xI M (,). (65) 

Comparing Eqs. (|55|) and (!65|) it is clear that, in the first segment, the B2-PLYP and Ai- 
B2-PLYP correlation integrands are formally identical. The only difference lies in the MP2 
correlation energies, which are not calculated with the same set of orbitals, as discussed pre- 
viously. A qualitative behaviour of the correlation integrand along the adiabatic connection 
is obtained when neglecting the variation of all terms that depend implicitly on v: 

rr-r-a x ,a c ,u -=(2)a x ,a c _ , ■. 

+2^ c LYP [^o]xX [Qxil] (z/). (66) 

A striking difference with the Ai-B2-PLYP correlation integrand (see Eq. (156]) ) is the positive 
slope in the second segment. This unphysical behaviour is directly related to our definition 
of the B2-PLYP integrand. It is a simple illustration of the fact that, when a c < a^, B2- 
PLYP does not rely on the density-functional perturbation theory we derived, by contrast 
to Ai-B2-PLYP. Still, after integration over the second segment, the B2-PLYP integrand 
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provides an energy contribution which differs from the Ai-B2-PLYP one as 

2(a x - V a x a cj ^ M p -^mp Jdz/ 
+ / X 2 (2a x - ^Jal - a c - 2i/)) £ C LYP [n* ] 

= 2(a x - a/^x - a c) \A X - a c x f^"*'" ~ ^S"* v ^^j , (67) 

if we neglect the variation of all terms that depend implicitly on v. As a result, B2-PLYP 
and Ai-B2-PLYP correlation energies will essentially differ by the MP2 term. 

Finally, we remark that since the B2-PLYP energy was not derived by consideration of 
the AC directly it is possible to construct a number of AC integrands for this approach. 
One alternative segmentation has already been presented in Eq. (I6ip . However, another 
possibility is to regard the B2-PLYP parameters as entirely empirical parameters, which 
simply scale the ACs derived for each component by a constant at all values of the interaction 
strength. A smooth AC integrand for B2-PLYP can then be obtained by summing these 
scaled components. However, whilst this integrand can be compared with the ab initio 
curves, the connection to the density-functional perturbation theory and the Ai-B2-PLYP 
methods presented here is lost. 




D. Summary 

A three-part segmentation of the exact exchange-correlation integrand has been proposed, 
which is directly connected to the double hybrid functionals of Ref. {40 1. Each segment of 
the AC has been expanded through second order within density-functional perturbation 
theory. When neglecting both second-order corrections to the density and density scaling, 
and using the Becke exchange functional in conjunction with the LYP correlation functional, 
the Ai-B2-PLYP integrand is obtained. An integrand expression has also been derived for 
the conventional B2-PLYP scheme. Both schemes are completely equivalent when a x = a c 
or, equivalently, Ai = A2. In this case, the second segment simply disappears. Interestingly 
for standard a x = 0.53 and a c = 0.27 values [yj], a x ~ 0.28 differs only by 0.01 from a c , as 



already pointed out by Sharkas et al. 39]. Still, since A 2 differs from Ai by a/o x — a c ~ 0.1 
(see Eq. (Jl8~]l ). the second segment represents 10% of the total AC which is not negligible. 
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TABLE I: The exchange-correlation integrands computed in this work ($° denotes the KS-BLYP 
determinant). 

Integrand Exchange Correlation Parameters 

ab initio Eq. ([M]) Eq. ([68]) 

Ai-B2-PLYP Eq. ([52]) Eq. §E§ Ai « 0.426 a , A 2 = 0.53 a 
Ax-DS-B2-PLYP Eq. ([52]) Eq. ([57]) Ai » 0.426, A 2 = 0.53 
B2-PLYP Eq. (JSU Eq. ([65]) a x = 0.53 fe , a c = 0.27 6 

BLYP Ar-1n |0 f 



a see Eq. 
b Ref. [lj 
c see Eq. (JED 

For comparison the methods used to determine accurate ab initio and pure density- 
functional estimates of the integrands have been outlined. All the schemes investigated in 
this work are summarized in Table [I] 



IV. RESULTS AND DISCUSSION 



In the present work we study the AC integrands for the species, H2, (He) 2, He-Ne, LiH, 
HF, N 2 and H2O. For H2 we consider the geometries R = 1.4 and 3.0 a.u., for the (He) 2 
and He-Ne van der Waals dimers we perform calculations at the equilibrium geometries 
of 5.612 a.u. 68] and 5.728 a.u. 69], respectively. For the remainin g fo ur molecular sys- 
tems LiH, HF, N 2 and H 2 we use the equilibrium geometries of Ref. [70] calculated at the 
CCSD(T)/cc-pVTZ level with all electrons correlated. All calculations of the AC integrands 
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73] using a modified version of the DAL- 



are performed in the aug-cc-pVTZ basis set 

TON2011 program [67], which contains implementations of the methodologies outlined in 
Section II I II and summarized in Table [B 
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A. ab initio ACs 



We have performed calculations using the methodology described in Section IIII Al 
for the species above. A range of interaction strengths in the interval v G [0,1], 
have been considered. In order to account for rapid curvature in the low-z/ part 
of the curve characteristic of statically correlated systems we have used the v values 
{0,10~ 6 , 10^ 5 ,10- 4 ,10- 3 ,10- 2 ,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1.0}. In Ref. [5l| a form 
for the AC integrand was proposed based on consideration of a simple two-state CI model 
and shown to have sufficient flexibility to reproduce correlation energies in systems exhibit- 
ing both static and dynamical correlation effects. Here we perform a least squares fit of this 
form 



W ,,AC-CI = _1+V5 ( 

4 



4(2 + y/5)a 2 + 5(3 + V5) 

(68J 



2^8(7 + 3y/E)a? + 16(2 + VE)asu + 10(3 + y/Z)s 2 v 2 

to the calculated ab initio estimates of [n] at the values of v above. The fitted values of 
the parameters a and s are reported in Table HT1 for each species in this study. Also reported 
is the quantity AE C = J* W^ AC ~ cl du - (Eqq - E nn - T s [n] - E ne [n] - E n [n] - E x [n]) which 
provides a consistency check for the quality of the fitted function as compared with the 
explicitly calculated correlation energy using non-interacting and interacting energies. As 
can be seen in Table [IT] these values are reasonable and of sufficient accuracy to allow these 
fitted functions to serve as a benchmark against which to compare double-hybrid integrands. 



B. The H2 molecule 



The H 2 molecule has been widely studied as a prototypical system, see e.g. Refs. [74j. |75| . 
exhibiting a smooth transition from predominantly dynamical correlation effects at short and 
equilibrium R values to predominantly static correlation effects a large R values. It has been 



argued by Gritsenko et al. 12j that for single hybrid exchange-based functionals the optimal 



fraction of orbital dependent exchange varies with R, approaching zero as R increases. More 
recently in Ref. j^li an analysis of the BLYP AC integrand , in comparison with the FCI AC 
integrand, showed that close to equilibrium R the BLYP functional provides a reasonable 
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TABLE II: Fitted coefficients a and s in Eq. (|68p for the species considered in this work. Also 
shown are the correlation energies calculated by integration of the fitted curves and the difference 
between these values and those calculated as AE C = J* Wc ' AC ~ CI df — (Eqq — E nn — T s [n] — E nc [n] — 
E u [n] - E x [n}) 



Molecule 


Q 


g 


Er 


AEr 


H 2 (R=1.4 a.u.) 


-0.171004 


-0.095425 


-0.039851 


6.91 x 10~ 6 


H 2 (R=3.0 a.u.) 


-0.153978 


-0.255931 


-0.076559 


-3.82 x 10~ 5 


(He) 2 


-0.513334 


-0.176830 


-0.079183 


5.62 x 10~ 6 


He-Ne 


-1.393157 


-0.806875 


-0.334652 


4.38 x 10~ 4 


HF 


-1.061605 


-0.775661 


-0.306331 


8.03 x 10~ 4 


LiH 


-0.220273 


-0.128740 


-0.053303 


6.37 x 10~ 5 


N 2 


-1.226312 


-1.201367 


-0.438569 


1.56 x 10~ 3 


H 2 


-0.993788 


-0.775566 


-0.301455 


8.37 x 10~ 4 



estimate of exchange and correlation energies. However for larger R values beyond ~ 7 
a.u. the estimate of the exchange energy is significantly too negative, whilst the correlation 
energy is significantly too positive. This is manifested by a much too flat shape for the 
BLYP correlation AC integrand at these geometries. At intermediate geometries R ~ 3 a.u. 
error cancellations between the exchange and correlation energies can lead to reasonable 
total energies. 

The exchange energy contributions to the two-parameter double hybrids are shown in 
Table ULTl The individual HF and density-functional type contributions are shown, calculated 
on the KS-BLYP $° determinant and n^o density, respectively. Their weighted contribution, 
of relevance to the double-hybrid functionals, is also tabulated. For comparison the exchange 
energies relevant to the ab initio estimates of the AC are also included, these are evaluated 
from the KS orbitals at v = which give the FCI density. For the H 2 molecule at R — 1.4 
a.u. it is clear that both the HF and density-functional estimates of the exchange energy are 
comparable. Their weighted average is also close to the exchange energy calculated for the 
KS orbitals giving the FCI density. At the longer bond length of R = 3.0 a.u. the difference 
between the HF and density-functional exchange contributions is much more pronounced. 
Comparing with the accurate FCI value in the same basis set it is clear that the density- 
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TABLE III: Exchange energy contributions for the Ai-B2-PLYP and B2-PLYP double hybrid 
schemes. The reference determinant $° is the KS-BLYP one. The HF exchange weight is set to 
a x = 0.53. For comparison the exchange energies £^ F [<3? KS ] of KS determinants constrained to 
yield accurate ab initio densities have also been included (see text for details). All values are given 
in atomic units. 









+(1- 


a x )E*[n^ } 




H 2 {R = 


1.4 a.u.) 


-0.6566 


-0.6563 


-0.6565 


-0.6608 


H 2 (R = 


3.0 a.u.) 


-0.4720 


-0.5061 


-0.4880 


-0.4769 


He 2 




-2.0295 


-2.0364 


-2.0327 


-2.0460 


HeNe 




-13.0517 


-13.1084 


-13.0783 


-13.0861 


LiH 




-2.1297 


-2.1292 


-2.1294 


-2.1369 


HF 




-10.3709 


-10.4404 


-10.4036 


-10.3870 


N 2 




-13.0855 


-13.1977 


-13.1382 


-13.0888 


H 2 




-8.9062 


-8.9674 


-8.9350 


-8.9149 



functional gives a much too negative exchange energy, as was also noted in Ref. [51J- Here 
we see that the weighted average used in the double-hybrid approaches significantly reduces 
this error. 

In Fig. [1] we present the correlation integrands for the double-hybrid approximations, 
as well as the BLYP and ab initio estimates. In the left-hand panel the correlation AC 



integrands for the methods considered are shown at R = 1.4 a.u. As was shown in Ref. [51 ] 
at this geometry the BLYP AC integrand is reasonable, though it tends to be too positive 
for larger v values. The challenge for double-hybrid approaches is to provide a model AC 
integrand which improves over the pure DFT integrand (in this case BLYP) whilst utilizing 
the DFT integrand where it is accurate. For the R — 1.4 a.u. geometry the total correlation 
energies in Table IIVI are all quite similar and close to the FCI estimate. This is also clear 
graphically from Fig. [TJ 

A number of differences between the methods become apparent when examining the 
correlation integrand models graphically. For B2-PLYP we see that the integrand in each 
segment is linear in the interaction strength. For the first segment this is because the 
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TABLE IV: The correlation energies for each segment of the double-hybrid ACs. For comparison 
the values for the correlation energies calculated from BLYP as well as accurate ab initio estimates 
of the AC are included. 
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FIG. 1: ACs for the H2 molecule at R = 1.4 a.u. (left panel) and R = 3.0 a.u. (right panel) calcu- 
lated using the conventional B2-PLYP functional (solid red line), the Ai-B2-PLYP two parameter 
double hybrid (green line), the Ai-DS-B2-PLYP functional (blue line) (which includes coordinate 
scaling contributions) and the standard BLYP functional (purple line). Also included is an accu- 
rate ah initio AC calculated at the FCI level (blue points). The function in Eq. (|68p has been 
fitted to this data (red dashed line) and the coefficients for this form are reported in Table HH 



integrand is dominated by a PT2 contribution based on a fixed set of orbitals. For the third 
segment this is because the uniform coordinate scaling contributions to the DFT correlation 
component are neglected. The most striking feature however is the intermediate interval 
where the B2-PLYP integrand is linear with positive slope and discontinuous at both Ai and 
A2. The significance of this section and the ambiguity in the choice of B2-PLYP AC were 
discussed in Sec. IIII CI The Ai-B2-PLYP variant also has an integrand consisting of three 
linear segments, however, continuity is restored at Ai, although a derivative discontinuity 
remains and a discontinuity is still present at A2. Note that by definition the Ai-B2-PLYP 
and B2-PLYP integrands are identical in the third segment. The crossing in the middle of 
the second segment, which implies that the Ai-B2-PLYP and B2-PLYP correlation energies 
obtained by integration are very close (as confirmed in Table IPVT ) . is consistent with the small 
difference between the AC lines in the first segment (see Sec. IIII CI) . The Ai-DS-B2-PLYP 
variant includes the effects due to uniform coordinate scaling. This affects the slope in the 
second segment and makes the integrand coincide with the BLYP integrand in the third 
segment. Interestingly, taking into account density scaling in the LYP correlation functional 
does not improve, in this particular case, the correlation energy when compared to FCI (see 
also Table [TV|) . 

In the right hand panel of Fig. [1] the same integrands are presented for H2 at R = 3.0 
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a.u. Here the behaviour of the integrands within each segment remains qualitatively similar, 
however, all of the models now give a significantly too positive correlation energy as can be 
seen from Table HVl The discontinuous behaviour at A2 also becomes much more pronounced. 
This can be understood by noting that (as shown previously in Ref. 51] ) the BLYP integrand 
is especially poor as the bond is stretched and static correlation becomes more important. 
It is clear therefore that any double-hybrid wishing to perform well for systems in which 
static correlation plays a significant role must be constructed in a different manner, no 
partitioning of the AC involving a pure density functional component would seem to be 
advantageous. Furthermore, the neglect of density scaling amounts to a linear approximation 
of the integrand and so may affect the accuracy of the model even in regimes where dynamical 
correlation is dominant. Note also, in the first segment, the difference between the Aj-B2- 
PLYP and B2-PLYP slopes which originates from the fact that the corresponding MP2 
correlation energies are calculated with different orbitals. These are mainly characterized by 
the fraction of HF exchange which is larger for B2-PLYP (0.53) than for Ai-B2-PLYP (0.426). 
The difference becomes substantial upon bond stretching. As expected from Sec. IIII CI the 
correlation energy in the second segment is then larger (in absolute value) for Ai-B2-PLYP 
than for B2-PLYP. This is graphically illustrated in the right panel of Fig. [TJ where the 
corresponding AC lines do not cross in the middle of the segment, like in the equilibrium 
geometry (left panel). 

Since the double hybrids considered are based on PT2 theories one of course should 
not expect their range of applicability to extend to strongly correlated systems. Nonetheless 
these considerations may be helpful if one wishes to design a double- hybrid functional capable 
of describing molecules over a reasonable part of their potential energy surfaces around the 
equilibrium structure. Furthermore, it highlights the point that for overall improvement of 
double-hybrid approaches one should carefully consider not only the nature of the density 
functional approximations employed but also the wave-function contributions. 

As is discussed in Appendix [B] we have chosen in the present work to evaluate the PT2 
contributions for the approximate functionals on a fixed set of orbitals. As a consequence 
the corresponding contributions to the AC are linear in the interaction strength. If the or- 
bital relaxation is taken into account at each interaction strength then the same integrated 
value of the energy would be obtained but the PT2 integrand would become curved due to 
a dependence of the orbitals on the interaction strength (see Ref for farther discussion). 
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To obtain higher accuracy from the wave-function contribution to the double-hybrid func- 
tional it would be desirable to introduce terms of higher order in the interaction strength 
without introducing significant extra computational cost. In this respect it is interesting 
to consider alternatives to PT2. Natural choices here would be higher order perturbation 
approaches or coupled-cluster type methodologies. However, the computational cost of these 
approaches is sufficiently high as to make them undesirable for application in this context. 
One interesting set of alternatives, which can be evaluated at a cost similar to that of PT2 
tey and , a S d lsCUSS e d by Furche fl do con taln higher order co—oos in V are the 
RPA correlation energies. Investigation of these variants of the correlation energy in the 
context of double-hybrid approaches based on a linear AC may be worthwhile. A number of 
empirical and range-separated approaches to combine DFT and RPA have already appeared 



in the literature, see for example, Refs. 
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82). 



C. (He) 2 and He-Ne van der Waals dimers 

In Fig. [2] we present the correlation integrands for two van der Waals dimers at their 
equilibrium geometries. The (He) 2 dimer has been widely studied as a prototypical system 



for examining van der Waals and dispersion interactions in DFT, see e.g. Refs. [52|, |83| 
and references therein. Methods which mix DFT with PT2 theory in a range-separated 
manner based on non-linear ACs have proven useful for the treatment of van der Waals 
and dispersion interaction energies. That the range separation of these interaction energies 
can be successful has recently been demonstrated by calculating ab initio estimates of the 
AC integrands along non-linear paths. However, for conventional double hybrid approaches 
such as B2-PLYP the description of dispersion interactions is far less satisfactory. Indeed 



empirical dispersion corrections have been developed to add to this functional 84J, despite 
its PT2 component. 

The exchange energy contributions for the (He) 2 and He-Ne van der Waals dimer systems 
are shown in Table HTT1 For the (He)2 dimer the HF and density-functional estimates of the 
exchange energy are similar and slightly more positive than the estimate based on KS orbitals 
giving the CCSD(T) density, which we will denote KS[CCSD(T)]. The weighted average of 
the exchange energy relevant to the double hybrids is therefore also reasonable. For the He- 
Ne dimer the HF estimate of the exchange energy is more positive than the KS[CCSD(T)] 
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FIG. 2: ACs for the van der Waals dimers (He)2 (upper left panel) and He-Ne (upper right 
panel) calculated using the conventional B2-PLYP functional (solid red line), the Ai-B2-PLYP two 
parameter double hybrid (green line), the Ai-DS-B2-PLYP functional (blue line) (which includes 
coordinate scaling contributions) and the standard BLYP functional (purple line). Also included 
is an accurate ab initio estimate of the AC calculated at the CCSD(T) level (blue points). The 
function in Eq. (|68|) has been fitted to this data (red dashed line) and the coefficients for this form 
are reported in Table HT1 The lower panels ((He) 2 left and HeNe right) show the interaction ACs 
for each method as defined in Eq. (f69|) . 

estimate, whilst the density-functional estimate is more negative. The weighted average is 
therefore much closer to the accurate value. 

The upper two panels in Fig.[2]show the total correlation ACs for (He) 2 and He-Ne respec- 
tively. The shape of these curves are similar in many respects to those observed for H2 near 
its equilibrium geometry. Although the standard BLYP functional now gives a too negative 
integrand at all interaction strengths. The general similarity between the (He)2, He-Ne and 
H2 (R — 1.4 a.u.) curves reflects the fact that on-atom dynamical correlation dominates the 
overall correlation energy contribution. Still, by contrast with H2, density scaling improves 
the correlation energy of both van der Waals dimers when compared to CCSD(T). These 
graphical results echo the observation by Sharkas et al. [39] when computing atomization 
energies with various double hybrid density-functionals; introducing density-scaling effects 
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does not systematically provide more accurate results when it is applied to approximate cor- 
relation density-functionals. This clear in the present work when comparing the left hand 
panel of Fig. [TJ and and upper left panel of Fig. [2J 

To examine the important correlation energy contributions to the interaction energies of 
the van der Waals dimers we have calculated the interaction ACs as in Ref. |34j. These are 
defined as 

W^friDimer, ™Atom 1, ^Atorn 2] = W£ [n D imer] ~ W^Atom l] ~ W c " [n At om 2], (69) 

where each of the atomic contributions are evaluated in the presence of the basis functions 
of the other atom, thereby accounting for the basis-set superposition error in the calculated 
difference. These integrands are presented for both systems in the lower two panels of Fig. |2j 
The trends in both figures are remarkably similar. The ab initio estimates of the interaction 
ACs show integrands which become smoothly more ne gat ive with increasing interaction 
strength. This is similar to the behaviour shown in Ref. [34| for the (He) 2 system at larger 
internuclear separations. The behaviour of the BLYP integrand is striking because it is 
significantly too negative in the low-z/ regime before switching to positive values at larger 
v values. This behaviour is also similar to that observed in Ref. 34J; that the curves for 
both systems are so qualitatively similar reflects the fact that their shape is determined 
mainly by their behaviour under uniform coordinate scaling according to Eq. l ISTl) . rather 
than the density on which they are evaluated. The interaction ACs for the double hybrids 
are significantly afflicted by the errors in the LYP contributions. In the first segment of the 
AC the influence of the density-functional component is clear, particularly for the density 
scaled variant, even though it is not the dominant contribution to the overall correlation 
energy AC. It is also notable that all of the double hybrid interaction ACs are significantly 
too flat in this section. In this respect, the double hybrid schemes do not seem to improve 
BLYP at all, with respect to the correlation energy, since this component becomes less 
attractive. Moreover, the exchange interaction energy which is repulsive at the BLYP level 
(+26 [iEh for (He)2 and +57 fiE^ for He-Ne) becomes attractive at the double hybrid level 
(-114 fiEh for (He) 2 and -146 fiEh for He-Ne) whilst the corresponding ab initio values (+74 
fiEh for (He)2 and +68 [lE^ for He-Ne) are clearly repulsive. Let us keep in mind though 
that along our approximate double hybrid AC (i) the reference density is the KS-BLYP one 
(which is therefore not affected by the MP2 treatment) and (ii) the density is not strictly 
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the same along the AC. Since the interaction energies considered here are very small, the 
density constraint might be important and contribute, through the orbital relaxation, to both 
exchange and correlation interaction energies. The density obtained at the double hybrid 
level should then be used as reference for setting up a true AC where the density constraint 
is indeed fulfilled. This should clearly be analyzed further in the future. Nevertheless, when 
comparing the total interaction energies, Ai-B2-PLYP is less repulsive (+29 \xEh for (He)2 
and +13.5 /iE h for He-Ne) than BLYP (+136.5 /iE h and +159 fiE h for He-Ne). The overall 
result is that the double hybrids based on a linear AC do too little to improve the description 
of both (He) 2 and He-Ne dimers. This rationalizes to some extent the need for empirical 
dispersion corrections even when using a functional such as B2-PLYP, see e.g. Ref. [84 1. 



D. LiH, HF, N 2 and H 2 molecules 

Finally in this section we examine a number of small molecular systems close to their 
equilibrium geometries. The exchange energy contributions for these systems are shown in 
Table IIHI For LiH both the HF and density-functional contributions are close to each other 
and reasonable compared to the KS[CCSD(T)] estimates. As a consequence the average 
relevant to the double-hybrid approximations is also of similar accuracy. For the HF, N 2 
and H 2 molecules the HF estimates of the exchange energies are reasonable, whilst the 
density- functional estimates are too negative in comparison with the KS[CCSD(T)] values. 
The averaging used in the double hybrid approaches therefore improves the exchange energy 
estimate relative to the pure density-functional approach. 

The correlation AC integrands for LiH, HF, N 2 and H 2 are shown in the four panels 
of Fig. |3j The plots for HF, N 2 and H 2 are qualitatively similar to those for the van 
der Waals dimers and the H 2 molecule at R = 1.4 a.u. The ab initio estimates of the 
correlation AC integrand are relatively subtly curved for each of these species, reflecting the 
shape expected for systems dominated by dynamical correlation close to their equilibrium 
geometries. The double-hybrid model correlation ACs for these systems follow similar trends 
to those discussed in Sections IIVBI and IIV CI It is perhaps noteworthy that the initial part 
of the ACs up to A x for these systems is reasonably well described by PT2 theory on the 
Ai-interacting system. 

The LiH molecule is more challenging for the DFT and double-hybrid models. The ab ini- 
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FIG. 3: ACs for the molecules LiH (top left panel), HF (top right panel), N2 (bottom left panel) 
and H2O (bottom right panel) calculated using the conventional B2-PLYP functional (solid red 
line), the Ai-B2-PLYP two parameter double hybrid (green line), the Ai-DS-B2-PLYP functional 
(blue line) (which includes coordinate scaling contributions) and the standard BLYP functional 
(purple line). Also included is an accurate ab initio estimate of the AC calculated at the CCSD(T) 
level (blue points). The function in Eq. (J68]) has been fitted to this data (red dashed line) and the 
coefficients for this form are reported in Table HH 

tio estimate is typical of a dynamically correlated system close to its equilibrium geometry. 
However, the BLYP functional gives a much too negative correlation integrand with a much 
too steep initial slope. This error is typical for many density functional approximations, 
which struggle to describe diatomic molecules composed of group 1 elements. These species 
typically have rather long equilibrium bond lengths compared to diatomic molecules com- 
posed of main group elements. In the first part of the model double-hybrid ACs it is clear 
that the Ai-interacting PT2 estimate is reasonable and is close to the ab initio estimate. In 
the intermediate region as more DFT contributions are included the correlation ACs for the 
Ai variants begin to slightly overestimate the correlation integrand. In the third section the 
double hybrid models inherit the large errors present in the density functional description of 
the correlation integrand for this species. The striking difference in behaviour between the 
correlation AC models for the LiH molecule and the other species considered here highlights 
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the difficulty in developing a truly transferable double-hybrid approach suitable for a wide 
range of systems, even when close to their respective equilibrium geometries. However, for 
a large number of species close to their equilibrium structures double-hybrid approaches are 
expected to be accurate. In future we hope that the AC analysis presented here can be 
extended to a much larger set of molecules and used to effectively evaluate double hybrid 
approaches based on a variety of wave-function and density-functional components. This 
approach could be effective in identifying and avoiding models which rely on large error 
cancellations and provide more stringent tests of the models than just post construction 
benchmarking against experimental data. 

V. CONCLUSIONS 

In this work we have explicitly derived the AC integrands underlying double-hybrid ap- 
proaches. The integrands have been calculated for the conventional B2-PLYP scheme and 
its so-called Ai variant which was obtained from second-order density-functional perturba- 
tion theory. These integrands were then compared graphically with benchmark ab initio 
estimates to assess their accuracy and a number of interesting features have been high- 
lighted. The approximate one- and two-parameter double hybrid ACs were divided into 
three segments, the first up to Ai being dominated by wave function theory contributions 
at the PT2 level, the second between Ai and A 2 involving both wave function and density- 
functional contributions and the final section beyond A2 involving purely density-functional 
contributions. 

Within each section of the approximate ACs the impacts of the approximations utilized 
in their derivations have been highlighted. In the first section the use of PT2 theory based 
on fixed orbitals gives a linear approximation to the AC and its slope can be understood 
from the nature of the orbitals determined for the Ai interacting system. In the intermediate 
region the behaviour of recently developed two-parameter forms sharply contrasts that of 
empirical forms, the latter giving ACs with positive slopes. In the final section, which is 
determined by density-functional contributions, the neglect of uniform coordinate scaling 
effects has been highlighted and leads to a linear approximation of the correlation AC. 
Inclusion of these effects can restore some curvature in the integrand, although the outright 
accuracy still depends heavily on the density-functional form employed. 
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The most striking feature of the approximate double hybrid correlation integrands is the 
presence of (derivative) discontinuities at the connecting Ai and A2 interaction strengths. 
Whilst these discontinuities do not affect the calculation of correlation energies (only right 



continuity is required 



51]). they may have implications for the determination and uniqueness 



of the multiplicative potentials associated with keeping the density constant along the AC. In 
future the design of double-hybrid models which avoid these features in their AC integrands 
should be considered. 

In the context of the van der Waals dimers the difficulties in constructing a double- 
hybrid approach based on the linear AC have been rationalized in terms of the interaction 
ACs. Here the failure of Ai-B2-PLYP and B2-PLYP approaches to account for longer- 
ranged interactions is evident and results in a significant underestimation of the correlation 
interaction energy. Recently, it has been shown that range-dependent generalized ACs can 
leverage physical insight about the range of interactions in this context to more effectively 
divide labour between the density- functional and wave function components 34|. We also 
note that the techniques employed in this work may be directly carried over to the analysis of 
range-separated double-hybrid methods by choosing an alternative non-linear AC integration 
path. 

Finally, we remark that the integrands presented here highlight that for more successful 
double-hybrid approaches it is essential to seek both improved wave function and density- 
functional components. In Section IIVBI we highlighted the RPA based methods as one 
possible route to include terms of higher-order in the interaction strength. Clearly it is a 
challenging task to introduce such higher-order contributions without incurring significantly 
increased computational cost. An equally challenging task is the construction of density 
functional components more compatible with these orbital based methodologies. It remains 
to be investigated if forms based on correlation functionals other than LYP can be more 
effective in this sense. The techniques used here could also be extended to further segment 
the AC in order to rationalize the behaviour of double hybrids with three or more parameters 
(see, for example, Refs. 20M22| ) . The ab initio estimates of the AC integrands may provide 
useful guidance in the development of these approaches and we expect that the type of 
analysis outlined here can play a central role in the future development of more robust 
double-hybrid approximations. 
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Appendix A: Integration of the exact segmented integrand 

Let us introduce the A-dependent decomposition of the exact ground-state energy 

Jx df 

= (V x \f + XW CC + V\V X ) + E^M, (Al) 
where the complement A-interacting Hxc density-functional energy equals 



EmM = / (*1Wee|^> dv 



A 

1 p\ 

(WlWnlW) du - (^|^ee|^) du, (A2) 
Jo 

this leads, according to Eqs. (jBJ) and OH]), to the expression given in Eq. (JIB"]) . In the exact 
theory, the energy E x does not depend on A. Still, it is convenient to derive its derivative 
with respect to A. It is obtained from the variational expression of the energy 

E x = mm {(# |T + XW CC + V\t>) + #Lm}> ( A 3) 
and the Hellmann-Feynman theorem which leads to 

^ = (* x \w ee \* x ) + ^M 

dA oX 

= (^ A |IU ee |^ A ) - E Hx [n*x] - A x [n^]. (A4) 

When integrating the segmented exchange-correlation integrand in Eqs. ( jl~Tj) and ( fl2l) . we 
therefore obtain from Eqs. (El), (GO), and (1A4I) . 

f Xl <\F, V 
E = E°+ ^du 
Jo du 

+ jf 2 U^IWeel^^-^Hxf^l-A^^lJ dlS, (A5) 



37 



which leads to Eq. ( !T4|) . 



Appendix B: Perturbation expansion of the exact segmented integrand 

According to Eqs. f lT5|) and (132]) . in the a = 1 limit, the energy E a,u ' v reduces through 
second order to 

E [2]u = (pu\f + u] y ee + V\&) + E^Jn^] + v 2 E® v . (Bl) 

From Eq. (TT8]) and the Hellmann-Feynman theorem, we obtain the first-order derivative 
expression 

= (<^\W ee \^) - ~ KM + —(v 2 E^y (B2) 

In addition, in the a = 1 limit, the first-order derivative in the third term on the right hand 
side of Eq. ( 133]) reduces through second order to 

J T?[2]\l,V 

—j— = (* A 'Hi/„|* J '>-£ Hx [n 1 ,» 1 ]-A^K, 1 ]+2A 1 £<S Al 

-/ dr (^i"-i + S |n *" i ) fa(2,A ' (r) ' (B3) 

Note that, in order to compute the MP2 term in Eq. (1B2|) . one would in principle need the 
response of the orbitals and their energies to the variations of v. Instead, we replace the 
//—dependent MP2 correlation energy by its value at u — \\ which, after integration over 
[0,Ai], gives the same result: 

A [ v 2 E$A du = 2uE^ du. (B4) 
o ^\ / Jo 

Combining all equations with Eq. ( ITT]) leads to the second-order expansion of the correlation 
integrand given in Eq. fl35l) . 

Appendix C: Correlation integrand associated to B2-PLYP 

Since E u ' u = E u ' u , the second and third term in the right-hand side of Eq. (162]) are derived 
exactly like in the Ai-B2-PLYP scheme, using the Hellmann-Feynman theorem. Similarly, 
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we obtain 

+ ^||(4-(a x -^) 2 )^ X ' a " (a ^ )2 )- (CI) 

In order to avoid the calculation of the orbital response to variations of u, we gather all MP2 
contributions as follows, 



o 



d / 2 ^ 

— I V /: ,„-, )<!// 

+ / , T- ^ E mp - K - K - v) )E Mp df 



rcix-y/ a x — a c 

which finally leads to the correlation integrand expression in Eq. (|65|) . 



f* ^ ^ 2uE^ ac du + 2(a x - v^r^)^ di/, (C2) 
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